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Рассмотрены два алгоритма конечноэлементного моделирования неоднородных свойств, в том числе при 
связности механических и электрических полей. В первом случае неоднородные механические и пьезоэлек- 
трические свойства тела задаются в виде аналитических зависимостей в некоторых классах Функции. Во 
втором предполагается, что Функции, описывающие неоднородные свойства, известны в дискретном наборе 
точек. Неоднородные свойства учитываются при вычислении элементных матриц. Для этого во втором слу- 
чае интерполяция неоднородных свойств внутри тела осуществляется двумя способами: с помощью сплай- 
нов, реализующих модель тонкой упругой пластины, и сплайнов Ж, Менге. Предложенные алгоритмы реали- 
зованы в виде модулей конечноэлементного комплекса АСЕГАМ. Разработанное программное обеспечение 
может быть использовано при решении статических и динамических задач для неоднородных упругих и 
электроупругих тел. Приведены соответствующие примеры. 
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Введение. Функционально-градиентные (ФГМ) материалы широко применяются в современном 
приборостроении, авиастроении, медицине и т. п. Физическое моделирование таких материалов 
дорого и не всегда возможно. Поэтому проектирование технических устройств, в которых исполь- 
зуются ФГМ, опирается на теоретическое исследование напряжённо-деформированного состоя- 
ния. В большинстве случаев оно основывается на численном математическом моделировании — 
например, в рамках метода конечных элементов (МКЭ). Поэтому разработка программного обес- 
печения (ПО), позволяющего проводить такое моделирование, представляется весьма актуаль- 
ной. В настоящей статье рассматривается создание такого ПО для конечноэлементного комплекса 
АСЕГАМ. 

В современных исследованиях выделяются два класса функционально неоднородных ма- 
териалов. В первом случае неоднородность возникает при поляризации керамики, когда элек- 
тродное покрытие имеет сложную геометрию [1]. В работе [2] рассматривается неоднородность 
другой природы. В неоднородном материале имеется трещина, которая подвергается структурно- 
му анализу при помощи метода Галёркина — Петрова. Вблизи трещины все свойства изменяются 
по одному закону (предполагается увеличение в Е?” раз). Исключение — тепловые коэффициен- 
ты, которые просто экспоненциально убывают по мере приближения к трещине. 

В настоящей статье разработан модуль конечноэлементного комплекса АСЕГАМ для мо- 
делирования функционально неоднородных материалов, в том числе пьезоэлектрических кера- 
мик с неоднородной поляризацией. При этом функция распределения неоднородных свойств 
материала, зависящая от положения точки в теле, предполагается известной в виде либо функ- 
циональной зависимости, либо дискретного набора значений величин соответствующих физи- 
ческих свойств в определённых узлах. Описана методика интерполяции этих свойств на исполь- 
зуемую конечноэлементную сетку. Приведены примеры решения задач для тел с неоднородны- 
ми свойствами. 


“Работа выполнена при финансовой поддержке РФФИ (гранты 13-01-00196-а, 13-01-00943-а, 13-01-12006-а). 
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Постановка задач для тел с неоднородными свойствами. Математическая модель состав- 
ного упругого, электроупругого и акустического тела с неоднородными свойствами твёрдых тел 
состоит из краевой задачи [3]: 

— для упругих и электроупругих тел 





р„@-а,р,и-У.“=Ё, У.09=0, (1) 
о=с". (Е+В,=)-е’ -Е, О-+с,О=е, -(=+6,2)+=,°.Е, (2) 
= = (Уи+ Уи") 2, Е=-Уф; (3) 
— для акустической среды 
р +У-М=0; м=Уф, р‚м=У- в; в=-М+ьЬум. (4) 


Р;с; 
Здесь в — тензор напряжений; р; — плотность тела; е — тензор деформаций; и — вектор пере- 
мещений; О — вектор электрической индукции; Е — вектор напряжённости электрического поля; 
Г, — вектор массовых сил; ф — электрический потенциал; аз, Ва» са — коэффициенты демпфиро- 
вания; с; , Е. =. — тензоры упругих констант, пьезомодулей и диэлектрических проницаемо- 
стей; индекс / обозначает номер тела в модели. В уравнениях для акустической среды: с — ско- 
рость звука; м — вектор скоростей; ц — потенциал скоростей; р — звуковое давление, Т — еди- 
ничный тензор. 

К уравнениям (1)—(4) добавляются механические, электрические и акустические гранич- 
ные условия. Тензоры с; , е, ь =., и плотность тел будем считать функциями от положения точ- 
ки в теле: 

рь =Ры(Х), < =* (х), =, ==, (х), 97 =е' (х). (5) 
Для учёта степени поляризации предлагается использовать следующую зависимость механиче- 
ских, электрических свойств материалов: 
— для тензоров с; и =, — 9=9'+Р(9°-9'); 


(6) 


— для тензора е, — ОР" 


Здесь через д обозначены компоненты соответствующих тензоров. При этом индексом 
Г обозначены компоненты тензоров для изотропного состояния, а индексом а — для анизотропно- 


го. Модуль вектора поляризации |Р] влияет на то, какие свойства будут проявляться в большей 
степени. Причём тензор пьезомодулей е’ будет нулевым для изотропных тел. 


Интерполяция неоднородных свойств. Реализация алгоритма учёта неоднородных свойств 
на этапе формирования конечноэлементных матриц предполагает, что при вычислении локаль- 
ных матриц известен набор коэффициентов, описывающих свойства материала. Предлагается два 
способа учёта неоднородности: свойства материала можно считать постоянными в пределах каж- 
дого конечного элемента либо вычислять набор необходимых коэффициентов для каждого узла 
разбиения. Для расчёта узловых значений коэффициентов разработан блок восстановления 
функций в узлах конечноэлементной сетки. 

Для блока восстановления функций реализованы несколько алгоритмов построения 
сплайнов. Во-первых, это &/7-р/ае зр/те [4], отвечающий модели бесконечно протяжённой пла- 
стины, деформируемой в некотором наборе точек изгибом. К решению задачи нас приведёт ми- 
нимизация функционала, описывающего полную свободную энергию изогнутой упругой пластины. 
Вид этого функционала: 


2 2 2 
+00 +00 0’ф 0°ф 0°ф | 
7, (Ф)= || Е 522% :| 59 ахау > тп. (7) 











—00 —© 
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Решение такой задачи ф(х‚‚у,)=Ё, /=1,№ имеет общий вид 


М 
ф( ху) = У си ПГ, + Син +СиоХ +Сы.зУ. (8) 
7=1 


Здесь 

2 2 

г =(х-х,) +(у-у,).. (9) 
Чтобы определить коэффициенты с; нужно решить следующую систему линейных алгебраических 
уравнений: 


м 

2 2 

Ф(хиу,) =>. пл, + См.1 + См+2Х; + См.зУ} =Р,, 
= 


= (хх, +(у,-У,), 1=1,М; (10) 


М М М 
а 
[= 7=1 


1=1 
Эта система не является положительно определённой (на главной диагонали — нули). Для её ре- 
шения используется метод Гаусса с выбором ведущего элемента. Плюс этого подхода — простота 
реализации. Однако существуют и более эффективные способы построения сплайнов. Рассмотрим 
теорию интерполяционных [”-сплайнов, которые позволяют получить функции с заданной степе- 
нью гладкости. Возникает следующая задача на минимизацию функционала гладкости: 











2, (Ф) > ти, © = феи" ("|9 (#') =, Г=ЬМ| , (11) 
1/2 
| 
где 7. (Ф)=|У”Ф а ГУ (29) ао! . (12) 


Здесь #', /=1,М — узлы интерполяции; Ё, /=1,М№М — заданные числа (значения восстанавлива- 


емой функции в узлах интерполяции). 

Вывод этого уравнения приводится в книге О. В. Ашкеназы [4]. При условии, что количе- 
ство точек интерполирования хотя бы вдвое превосходит заданную степень гладкости /7 (в рас- 
сматриваемом случае это выполняется практически всегда), существует единственное решение 
задачи (12), общий вид которого: 


ФЕ =УсЕ, (ЕЕ) +р(Е), ЕЕК". (13) 
[= 
С точностью до положительного множителя: 


о” и ни пи -и|, при чётном п, 


Е (и, и) = (п-1)/2 2т-п .. (14) 
(-1) ии", при нечётном п. 
Второе слагаемое является многочленом степени т-1: 
М 
р(Ё) => а,Ф, (Р)еР,у. (15) 
У=1 
Причём ф, (Е), У=1, М, М=С”, — это различные одночлены степени, меньшей или равной 
т. Чтобы определить коэффициенты сплайна, нужно решить следующую систему уравнений: 
Ус (Е) + У аф, (ЕЁ =’, Г=ЕМ; 
Е м=1 (16) 


Ус, (#7) =0, у= 1, М. 
7 
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Эта система решается с помощью алгоритма Ж. Менге. Некоторое подмножество точек выделяет- 
ся для построения системы фундаментальных многочленов. Её порядок зависит от выбранной 
степени гладкости сплайна. Таким образом исходная система разделяется на две: первая (мень- 
шего размера) — для нахождения системы многочленов и основная — для нахождения сплайна. 
Основная матрица в этом методе является положительно определённой и симметричной, что поз- 
воляет использовать быстрый и устойчивый метод квадратного корня для её решения. 
Программная реализация в АСЕГАМ. Модуль построения сплайнов — отдельный исполняемый 
файл с собственным интерфейсом, позволяющим использовать его для различных целей. Связь 
между частями комплекса осуществляется при помощи набора утилит. Это позволяет заменять 
отдельные части пакета независимо друг от друга. Для генерации входных данных, описывающих 
неоднородные свойства, была создана специальная графическая оболочка. В ней можно задать 
неоднородности для конкретной модели, построенной в АСЕГАМ, и визуализировать их. Кроме 
этого проводится проверка введённых пользователем данных на непротиворечивость физической 
модели материала. 

Чтобы пользователю было удобно оценивать поведение восстановленной с помощью 
сплайнов функции, разработан визуализатор на основе тепловых карт. К тепловой карте прила- 
гаются отметки о максимальном и минимальном значении функции — им соответствуют опреде- 
лённые цвета. Карта строится на основе конечноэлементной сетки, поэтому на границах элемен- 
тов могут возникать артефакты изображения. Это происходит при относительно большой разнице 
значений функции между двумя элементами. На этот случай предусмотрена возможность сглажи- 
вания изображения на основе классического фильтра Гаусса. Все представленные утилиты реали- 
зуют программные интерфейсы для взаимодействия с данными в форматах АСЕГАМ и снабжены 
инструментами оценки эффективности модулей — записывается время выполнения всех опера- 
ций, сложных с точки зрения вычислений. 

Пример визуализации неоднородности представлен на рис. 1. 


Значения 
У. 7АЦЕНОО] Е 


9.4103Е—001 и 


Полябизация | 


Сглалить 





Рис. 1. Визуализация неоднородных свойств 


Прямые задачи для неоднородных тел. Рассмотрим работу описанного ПО на примере вдав- 
ливания штампа (на верхнюю сторону которого приложено равномерное давление) в упругую не- 
однородную полосу. Размеры штампа — 0,4 х 0,1 см. К верхней границе приложены равномерно 
распределённые нормальные напряжения (интенсивность — 103 Н/м”). Материал штампа — изо- 
тропный модуль Юнга 2,0 х 10" Н/м”. Коэффициент Пуассона — 0,29. Размеры полосы — 
3х 0,5 см. Нижний край закреплён. Материал — цирконат-титанат свинца (Р77-4). Упругие свой- 
ства полосы изменяются по закону д = д ехр(29х,!), где д — компоненты тензора упругих посто- 
янных неоднородного тела, д — соответствующие компоненты для однородно поляризованной 
керамики. 
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На рис. 2 представлена геометрия задачи. Цветом выделяется неоднородность упругих 
свойств и распределение вертикального смещения на деформированном состоянии области. 


ХЗ 





3 Х! (см) 


ХЗ 





у 3 Х! (см) 


6) 


Рис. 2. Вдавливание штампа в неоднородную полосу: а — исходная геометрия задачи с тепловой картой неоднородной 
жёсткости, 6 — деформированное состояние пластины с распределением вертикальной компоненты смещения 


На рис. 3, Ч и 5 представлены расчёты по излучению акустических волн в воду (область 
серого цвета на рис. 3, изображена только часть) с помощью пьезоэлемента в виде полупассивно- 
го биморфа. Верхняя часть биморфа изготовлена из меди и имеет длину 0,125 м при толщине 
0,011 м. Нижняя часть состоит из керамики ЦТБС-3. Её длина 0,1 м, толщина — 0,01 м. Левый 
край биморфа закреплён. На верхней и нижней грани керамической пластинки заданы электриче- 
ские потенциалы, образующие разность в 100 В. 





Рис. 3. Визуализация неоднородности по электроупругим свойствам 


Рассматривался случай неоднородной поляризации пьезокерамической пластины, при ко- 
тором левый край поляризован полностью, а правый не поляризован (см. рис. 3). Модуль вектора 


поляризации |Р| в соотношениях (6) изменяется по закону |Р|(х) =-9х +0,1. 


На рис. 4 представлено распределение потенциала скоростей в жидкой среде при часто- 
те 10 кГц. 





а) 6) 


Рис. 4. Поле потенциала скорости в жидкости для неоднородной среды: а — однородный случай, 6 — неоднородный 
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На рис. 5 изображена амплитудно-частотная характеристика (АЧХ) вертикального смеще- 
ния крайнего левого верхнего узла пластины. 


и у——— ще 





10 1] 12 13 14 15 
Частота. кГц 
а) 


10 1] 12 13 14 15 
Частота, кГц 
6) 


Рис. 5. Результаты решения задачи для излучателя в акустической среде. 
АЧХ крайнего правого верхнего узла излучателя: а — для однородной среды, 6 — для неоднородной среды 


Заключение. Рассмотрены алгоритмы конечноэлементного моделирования неоднородных 
свойств упругих и электроупругих тел, а также их программная реализация в АСЕГАМ, Примене- 
ние предложенных моделей позволяет решать прямые задачи расчёта напряжённо-дефор- 
мированного состояния различных устройств и приборов. Кроме того, представленные методы и 
программы (например, в сочетании с генетическим алгоритмом) можно использовать при реше- 
нии задач оптимизации конструкций и обратных задач идентификации неоднородных свойств, 
возникших в результате технологических или эксплуатационных процессов. 
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ЕЕ 5ТМОГАТТОМ ОЕ ЗТКОСТОВЕ$ ЕКОМ МОМНОМОСЕМЕОЦ$ МАТЕКТАТ$ 
МЛТН СОМРЫСАТЕО РВОРЕВТТЕ$” 


Р. А. Одапежуам, А. М. $о!оууеу 
(Роп Чае Тесптса! Упмег®Ку) 


Тио Нпйе еетепЕ хти/авоп а/догт5 Гог Веегодепеоиз ргорегйез тсиатд пе соппеснуийу оЁ тесваги/са!/ апа 
е/есиса! Пе/а® аге сопз!аегец. [п {те Нг5Е пхапсе, Веегодепеоиз теспапгса! апа р/егое/есй"с ргорегйез ог! те 
роду аге дмеп т {пе Гогт о! апа/уйса! ге/аноп$/рз т себат с/аз5ез о! Гипсйоп$. Тйе 5есопа тетоа а5зитез 
паЁ те ипсйоп$ аезспбта пе Реегодепеои$ ргорегИе$ аге Кпоип т а @осгее $еЁ о! рот. Тпе свапдез$ 
саижеа Бу пвотодепеои$ ргорегИе$ аге {аКеп тю ассоипЕ ипаег те са/сшавоп о! {Ве еетепЕ тайтсе®. Рог 5, 
1 Пе 5есопа сае, Пе пщегро/айоп о! Веегодепеоиз ргорегйез иййт {те Воду 15 ипретепед т о игау: и5- 
па т-р/ае 5р/тес, апа Л. МетдиеЕ5 $р/тез. Тре ргорозеа а/догЁНт$ аге ипретепЕеа аз тодше$ о! {Ве АСЕ- 
[АМ Нппе еетепЕ сотрЕх. Ехатр/ез о! изта те деуеюореа 5оКугаге Гог {Пе сасш!айоп о! Те $е55 {ат $аЁе 
ог ппотодепеои$ ева5Ёс апа еесгое/азИс БосГе$ аге ргезегед. 

Кеуигога$: Гипсбопа/у тпвотодепеои$ е/азНс апа е/есгоев5йс тайепа/5, $р/те$, Ипе-еетепЕ тетпод. 


“ ТНе гезеагсн 15 Чопе мИЛ {Не Япапса! зиррог гот ВЕЕТ (дгап5 13-01-00196-а, 13-01-00943-а, 13-01-12006-а). 
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